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The  main  accomplishments  with  the  technical  details  have  been  reported  in  the  archival 
publications.  Therefore,  below  we  will  list  the  main  findings  with  the  corresponding  figures  from 
our  presentation  at  the  2015  Program  Review  Meeting  and  short  explanations.  We  also  include 
the  references  to  our  publications  for  the  derivations  and  the  detailed  explanations. 

Main  accomplishments 

1.  A  two-stage  time-integration  technique  for  elastodynamics  and  acoustic  wave 
propagation  problems 

The  application  of  the  finite  elements  or  any  other  space-discretization  method  to  the  partial 
differential  equations  of  linear  elastodynamics  and  acoustics  leads  to  a  system  of  ordinary 
differential  equations  in  time,  see  Fig.  1 . 

Semi-discrete  formulation  of  linear  elastodynamics 

Space  discretization  methods: 

a)  Finite  element  method 

b)  Spectral  elements  method 

c)  Isogeometric  elements 

d)  Boundary  element  method 

e)  smoothed  particle  hydrodynamics  (SPH)  method 

f)  others 

Semi-discrete  MU  +  CU  +  KU  =  R  (1) 
system  of  equations  t/(,  =  0)  =  f/05  =  0)  =  Ko 

1  The  exact  time  integration  of  the  semi-discrete  system,  Eq.  (1 ),  may  yield  very 
inaccurate  solutions  of  the  original  (before  the  space  discretization)  system  of 
elastodynamics  equations  due  to  large  spurious  high-frequency  oscillations. 

2.  The  goal  of  calculations  is  not  an  accurate  integration  of  Eq.  (1 )  but  accurate 
solutions  of  original  elastodynamics  problems  before  the  space  discretization. 

3.  Accurate  time  integration  of  low  modes  and  filtering  spurious  high  modes  are 
needed  ! 

Fig.  1 
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Issues  with  existing  time-integration  methods 
for  wave  propagation  problems  due  to 
spurious  high-frequency  oscillations 

□  For  wave  propagation  and  impact  problems,  existing  time-integration  methods  may 
yield  locally  divergent  results  at  small  time  increments  and  mesh  refinement  due  to 
large  spurious  high-frequency  oscillations. 

□  There  is  no  quantification  of  the  range  of  spurious  oscillations  for  semi-discrete 
elastodynamics  equations 

□  There  is  no  understanding  how  to  select  the  amount  of  numerical  dissipation  (or 
artificial  viscosity)  and  the  size  of  time  increments  with  existing  time-integration 
methods  for  wave  propagation  problems 

“In  no  way  can  today's  computer  program  for  wave  propagation  and  impact  be  treated  as  'black  boxes'.  A 
minimum  6  months  to  2  years  of  experience  is  needed  to  be  able  to  use  such  programs  successfully" 

(Concept  and  applications  of  finite  element  analysis,  R  D.  Cook  et  al,  2002). 

□  Usually,  accuracy  analysis  of  existing  time-integration  methods  with  numerical 
dissipation  (or  artificial  viscosity)  does  not  include  the  effect  of  the  number  of  time 
increments  on  the  accuracy  of  final  numerical  results. 

•  Numerical  dissipation  used  in  existing  time-integration  methods  for  filtering  the 
spurious  high  modes  affects  the  low  modes  of  a  numerical  solution  resulting  in 
inaccurate  solutions,  especially  for  long-term  integrations  (the  spectral  radius  r" 
should  be  close  to  1  for  actual  frequencies;  if  r  =  0.99  for  n  =  1 ,  then  rn  =  0.99100  = 

0.366  for  n  =  100,  or  if  r  =  0.999  for  n  =  1,  then  r"  =  0.9991000  =  0.3677  for  n  =  1000) 

□  The  selection  of  the  most  effective  time-integration  method  for  wave  propagation 
problems  among  known  methods 

Fig.  2 

The  main  issues  with  the  integration  of  this  system  of  ordinary  differential  equations  are  listed  in 
Fig.  2.  The  most  important  issue  is  related  to  divergent  results  at  time  and  mesh  refinements 
(e.g.,  below  see  Fig.  7  with  no  filtering).  To  resolve  these  issues  we  have  developed  a  two-stage 
time-integration  approach  consisting  of  the  stage  of  basic  computations  and  the  filtering  stage. 
The  idea  of  the  two-stage  time-integration  procedure  is  very  simple.  Because  for  the  ordinary 
differential  equations  of  linear  elastodynamics  and  acoustics  there  is  no  interaction  between 
different  modes  during  time  integration  (due  to  modal  decomposition  they  are  integrated 
independently  and  do  not  affect  one  another,  see  Fig. 3),  the  most  accurate  time-integration 
method  (without  numerical  dissipation  or  artificial  viscosity)  should  be  used  for  basic 
computations  (one  stage),  especially  for  a  long-term  integration.  This  means  that  all  modes 
(including  high  modes)  are  integrated  very  accurately,  and  the  solution  includes  high-frequency 
oscillations.  Then,  for  damping  out  high  spurious  (inaccurate)  modes,  a  time-integration  method 
with  large  numerical  dissipation  (or  with  artificial  viscosity)  can  be  used  as  a  pre-,  or  post¬ 
processor  (another  stage).  This  method  can  be  considered  a  filter  of  high  modes;  see  Fig.  4. 
Usually,  a  small  number  of  time  increments  (we  always  use  5  positive  plus  5  negative  time 
increments)  is  necessary  for  the  pre-,  or  post-processing  (filtering)  stage,  with  no  error 
accumulation  at  low  modes.  This  technique  yields  no  error  accumulation  due  to  numerical 
dissipation  (or  artificial  viscosity)  at  the  stage  of  basic  computations  and  does  not  require  any 
guess  for  the  selection  of  numerical  dissipation  or  artificial  viscosity  as  in  the  existing 
approaches.  One  of  the  impotent  components  of  the  new  numerical  approach  is  the  quantification 
of  the  range  of  spurious  frequencies.  Using  the  dimensionless  analysis  of  the  1-D  impact  wave 
propagation  problem  (all  frequencies  are  excited  for  this  problem),  defining  the  amplitudes  of 
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spurious  frequencies  as  the  difference  between  the  exact  and  numerical  velocities  before  and 
after  the  wave  front  for  the  1-D  problem,  we  have 

A  new  two-stage  time-integration  approach  for 
eiastodynamics  problems 
(basic  computations  and  the  filtering  stage) 

Semi-discrete  ^0  +  CLJ  +  KU  =  R 
system  of  equations  ^  =  0)  =  U0,  U(t  =  0)=V0 

Accurate  time  integration  of  low  modes  and  filtering  spurious  high  modes  are  needed  ! 

J3L  Modal  decomposition 

ii,  +  2  %-u,  +  co~ut  =rt,  i  =  1,2 

Hi=f(U0),  u°  =  f2(V0),  rt  =  f3(R) 

-U- 

All  modes  are  integrated  independently.  There  are  no 
interaction  between  modes  during  time  integration 

(spurious  modes  do  not  affect  other  modes) 

There  is  no  necessity  to  filter  spurious  oscillations  at  each  time  increment. 

Spurious  high-frequency  modes  can  be  damped  out  at  the  end  of  calculation  ! 

Fig.  3 


A  new  two-stage  time-integration  approach  for  eiastodynamics 
(basic  computations  and  the  filtering  stage) 

A  New  Two-stage  Time-integration  Procedure  ([1]) 

1.  The  basic  computations,  especially  for  a  long-term  integration,  can  be  done  by  a  high- 
order  time-integration  method  with  zero  (or  very  small)  numerical  dissipation  and  small 
sizes  of  time  increments.  This  means  that  ail  modes  are  integrated  very  accurately. 

The  trapezoidal  rule  is  best  method  for  basic  computations  among  all  implicit  second-order 

methods  III 


2.  Then,  for  the  damping  out  of  spurious  high  modes,  a  time-integration  method  with  large 

numerical  dissipation  can  be  used  for  a  small  number  of  time  increments  as  a  pre-  or  post¬ 
processor  (the  filtering  stage).  This  method  can  be  considered  as  a  filter  of  spurious  high 
modes. 

Two  different  time-integration  methods  are  used  with  the  two-stage  time-integration  approach  I 

Existing  Approaches 

A  time-integration  method  with  large  numerical  dissipation  at  high  modes  used  for  each  time 
increment,  leads  to  accumulation  of  a  numerical  error  for  low-mode  terms  as  well, 
especially  at  the  long-term  integration 


Fig.  4 
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found  the  analytical  expression  for  the  range  of  spurious  frequencies  in  terms  of  the  size  of  time 
increments  at  the  filtering  stage.  Because  the  range  of  spurious  frequencies  depends  on  the  mesh 
size  and  the  observation  time  and  is  independent  of  the  initial  and  boundary  conditions,  the 
analytical  expression  for  the  range  of  spurious  frequencies  is  applicable  to  all  elastodynamics 
and  acoustics  problems.  This  expression  includes  two  parameters  ai  and  ao  that  depend  on  the 
type  and  the  order  of  the  space  discretization  method;  see  Fig.  5.  The  values  of  the  parameters  ai 
and  a2  for  the  low-  and  high-order  finite  elements,  spectral  elements,  isogeometric  elements  and 
the  linear  finite  elements  with  reduced  dispersion  are  determined  and  shown  in  Fig.  6.  It  should 
be  mentioned  that  the  two-stage  time-integration  approach  with  the  filtering  stage  yields  accurate 
convergent  results  at  time  and  mesh  refinements;  e.g.,  see  the  numerical  results  for  the  2-D 
impact  problem  in  Fig.  7.  Without  the  filtering  stage,  the  error  in  velocities  (or  stresses)  may 
exceed  300%  and  more;  see  Fig.  7. 
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Coefficients  ay  and  a2  for  different 
space-discretization  methods  [1] 

(used  for  the  calculation  of  the  time  increments  at  the  filtering  stage) 


At  =  a(Nx  )AyQ(»j(A)  a(Nx )  = 


yAxj 


=  a, 


f  cT^“2 


Types  of  Elements 

Coeffi- 

dents 

I5' 

^nd 

3“ 

4™ 

5“ 

10“ 

Standard  Finite  Elements 

<*1 

0.3574 

0.3156 

0.4485 

0.5495 

0.5115 

(Consistent  mass  matrix) 

a  2 

0.3204 

0.2364 

0.1571 

0.1111 

0.1078 

Spectral  Elements 

al 

0.3342 

0.448 

0.5659 

0.4790 

0.4461 

0.4317 

(Lmnped  mass  matrix) 

a2 

0.3363 

0.1845 

0.1139 

0.1128 

0.1073 

0.0759 

Isogeometric  Elements 

0.2513 

0.2311 

(Consistent  mass  matrix) 

°2 

0.2035 

0.1508 

Linear  Elements  with  Reduced 

«1 

0.2979 

Dispersion  (Non-diagonal  matrix) 

a? 

0.2074 

Linear  Elements  with  Reduced 

Or 

0.3296 

Dispersion  (Diagonal  matrix) 

a?. 

0.2180 

[1]  A.  V  Idesman,  D.  Pham,  Accurate  solutions  of  wave  propagation  problems  under  impact  loading  by  the 
standard,  spectral  and  isogeometric  high-order  finite  elements.  Comparative  study  of  accuracy  of  different  space- 
discretization  techniques.,  Finite  Elements  in  Analysis  and  Design,  88,  2014,  pp.  67-89. 


Fig.  6 


2D  impact  axisymmetric  problem:  Convergence  at  mesh 
refinement  using  the  new  technique  (with  the  filtering  stage) 


B  C 


V0 


r 


Without  the  filtering  stage 


z  1 .  Mesh  40  X  100  =  4000  Q4  elements 

2.  Mesh  80  X  200  =  16000  Q4  elements 

3.  Mesh  160  X  400  =  64000  Q4  elements 

4.  Reference  Solution 


Velocity  distribution  along  axis  of  revolution  for  different  meshes  with  linear  elements 


Fig.  7 
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Our  publications  with  a  detailed  description  of  the  results  presented  above. 

1 .  Idesman  A.  V.  ,  Accurate  finite  element  modeling  of  wave  propagation  in  composite  and 
functionally  graded  materials,  Composite  Structures,  2014,  1 17,  pp.  298-308. 

2.  Idesman  A.V.,  Pham  D.,  Foley  J.,  Schmidt  M.  Accurate  solutions  of  wave  propagation 
problems  under  impact  loading  by  the  standard,  spectral  and  isogeometric  high-order 
finite  elements.  Comparative  study  of  accuracy  of  different  space-discretization 
techniques,  Finite  Elements  in  Analysis  and  Design,  2014,  88,  pp.  67-89. 

3.  Idesman  A.V.,  Pham  D.  ,  Accurate  finite  element  modeling  of  acoustic  waves,  Computer 
Physics  Communications,  2014,  185,  pp.  2034-2045. 

4.  Idesman  A.V.,  Mates  S.  P.,  Accurate  finite  element  simulation  and  experimental  study  of 
elastic  wave  propagation  in  a  long  cylinder  under  impact  loading,  International  Journal  of 
Impact  Engineering,  2014,  71,  pp.  1-16. 

5.  Idesman  A.V.,  Pham  D.,  Finite  element  modeling  of  linear  elastodynamics  problems 
with  explicit  time-integration  methods  and  linear  elements  with  the  reduced  dispersion 
error,  Computer  Methods  in  Applied  Mechanics  and  Engineering,  2014,  271  pp.  86-108. 


2.  New  finite  elements  with  the  reduced  dispersion  error  for  explicit  time-integration 
methods 

We  have  developed  two  finite  element  techniques  with  reduced  dispersion  for  linear 
elastodynamics  that  are  used  with  explicit  time-integration  methods.  These  techniques  are  based 
on  the  modified  integration  rule  for  the  mass  and  stiffness  matrices  and  on  the  averaged  mass 
matrix  approaches  (see  Fig.  8)  that  lead  to  the  numerical  dispersion  reduction  for  linear  finite 
elements  from  the  second  order  to  the  fourth  order.  The  analytical  study  of  numerical  dispersion 
for  the  new  techniques  has  been  carried  out  in  the  1-D,  2-D  and  3-D  cases.  The  numerical  study 
of  the  efficiency  of  the  dispersion  reduction  techniques  includes  two-stage  time-integration 
approach  with  the  filtering  stage  that  quantifies  and  removes  spurious  high-frequency 
oscillations  from  numerical  results.  We  have  found  that  in  contrast  to  the  standard  linear 
elements  with  explicit  time-integration  methods  and  the  lumped  mass  matrix,  the  finite  element 
techniques  with  reduced  dispersion  yield  more  accurate  results  at  small  time  increments  (smaller 
than  the  stability  limit)  in  the  2-D  and  3-D  cases.  The  recommendations  for  the  selection  of  the 
size  of  time  increments  have  been  suggested.  The  new  approaches  with  reduced  dispersion  can 
be  easily  implemented  into  existing  finite  element  codes  and  lead  to  significant  reduction  in 
computation  time  by  a  factor  of  10  -  1000  and  more  compared  with  the  standard  finite  element 
formulations  at  a  given  accuracy.  We  should  also  mention  that  the  fonnulations  with  reduced 
dispersion  and  the  explicit  time-integration  methods  can  be  directly  and  efficiently  used 
(practically,  without  modifications)  on  parallel  computers  because  the  numerical  algorithm 
includes  the  matrix  and  vector  multiplications  and  can  be  perfonned  at  the  element  level  without 
the  solution  of  a  system  of  algebraic  equations. 
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Modification  of  FE  equations  for  explicit  time- 
integration  methods 

MU  +  KU  =  R  is  replaced  by  DV  +  KU  =  R ,  MV  =  DU 

or 

dO  +  MD  'KU  =  MD  R 

From  the  dispersion  analysis  (the  central  difference  method): 

□  Averaged  mass  matrix  technique  ([1]) 

/  =  (3  -  r 2 )/  2 


□  Modified  integration  rule  technique  ([1]) 

aAI  =  V(4-r2)/3, 


ak 


4(2u-l) 

3(4v-3) 


where  z=c2kt/dx,  C2=J 


E 


2p(\  +  v) 


E,  V,  p  are  the  material  parameters 


[1]  Finite  element  modeling  of  linear  clastodynamics  problems  with  explicit  time-integration  methods  and 
linear  elements  with  the  reduced  dispersion  error,  Idesman  A.V.,  Pham  D.,  Computer  Methods  in  Applied 
Mechanics  and  Engineering.  271, 2014.  pp.  86-108. 


Fig.  8 


Our  publications  with  a  detailed  description  of  the  results  presented  above. 


1.  Idesman  A.V.,  Pham  D.,  Finite  element  modeling  of  linear  elastodynamics  problems 
with  explicit  time-integration  methods  and  linear  elements  with  the  reduced  dispersion 
error,  Computer  Methods  in  Applied  Mechanics  and  Engineering,  2014,  271  pp.  86-108. 

2.  Idesman  A.V.,  Pham  D.,  Foley  J.,  Schmidt  M.  Accurate  solutions  of  wave  propagation 
problems  under  impact  loading  by  the  standard,  spectral  and  isogeometric  high-order 
finite  elements.  Comparative  study  of  accuracy  of  different  space-discretization 
techniques,  Finite  Elements  in  Analysis  and  Design,  2014,  88,  pp.  67-89. 

3.  Idesman  A.V.,  Pham  D.  ,  Accurate  finite  element  modeling  of  acoustic  waves,  Computer 
Physics  Communications,  2014,  185,  pp.  2034-2045. 


3.  Comparison  of  accuracy  of  different  space-discretization  techniques  used  for  wave 
propagation 

One  issue  with  numerical  solutions  of  wave  propagation  problems  solved  by  different  space- 
discretization  techniques  is  the  presence  of  spurious  high-frequency  oscillations.  The  range  of 
spurious  high  frequencies  is  different  for  different  space-discretization  methods.  We  have  not 
seen  in  the  literature  the  numerical  approaches  that  quantify  and  filter  out  all  spurious 
oscillations  from  numerical  solutions  even  solved  by  the  popular  space-discretization  techniques 
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such  as  the  standard,  spectral  and  isogeometric  high-order  finite  elements.  We  have  resolved  this 
issue  by  the  application  of  the  two-stage  time-integration  approach. 

Another  important  result  is  related  to  the  comparison  of  accuracy  of  different  space- 
discretization  techniques  used  for  the  solution  of  transient  acoustics  and  elastodynamics 
problems.  One  way  to  compare  the  accuracy  of  these  techniques  is  based  on  the  evaluation  of  the 
numerical  dispersion  error.  However,  this  error  is  defined  for  one  selected  mode  of  the 
elastodynamics  system  of  equations  and  does  not  prescribe  the  combined  effect  of  all  modes  on 
the  accuracy.  Moreover,  the  dispersion  error  does  not  estimate  the  effect  of  the  observation  time 
on  the  accuracy  of  numerical  results.  Therefore,  for  the  comparison  of  the  accuracy  of  different 
techniques,  we  solved  the  1-D  impact  problem  for  which  all  modes  are  excited  and  which  has  a 
very  simple  analytical  solution  at  short  and  long  observation  times.  We  should  also  mention  that 
because  the  computational  costs  of  different  space-discretization  techniques  at  the  same  number 
of  degrees  of  freedom  (dof)  are  different,  then  it  is  necessary  to  compare  the  efficiency  of  these 
techniques  by  the  estimation  of  the  computational  costs  at  a  given  accuracy.  The  findings  can  be 
summarized  as  follows. 

The  two-stage  time-integration  technique  yields  accurate  numerical  results  for 
elastodynamics  problems  solved  with  different  space-discretization  approaches.  For 
example,  due  to  the  spurious  oscillations,  we  have  not  seen  in  the  literature  the  accurate 
numerical  solutions  of  elastodynamics  problems  at  impact  loading  (or  high-frequency 
loading)  for  the  spectral  and  isogeometric  elements  as  well  as  for  the  standard  high-order 
finite  elements  even  in  the  1-D  case.  The  applicability  of  the  new  approach  to  different 
space-discretization  techniques  is  based  on  the  fact  that  for  all  these  techniques,  lower 
frequencies  are  resolved  more  accurately  than  higher  frequencies.  Then,  by  the 
quantification  of  the  inaccurate  high  frequencies  (in  terms  of  the  coefficients  ai  and  ?lt) 
and  by  their  filtering  at  the  filtering  stage,  we  obtain  accurate  numerical  results  without 
the  spurious  oscillations;  see  Figs.  9-10  (after  filtering).  The  new  coefficients  &\  and  a.2 
define  different  ranges  of  the  spurious  (inaccurate)  frequencies  for  different  space- 
discretization  methods  and  are  used  for  the  filtering  of  the  spurious  oscillations  in 
numerical  solutions  of  the  1-D  as  well  as  multi-dimensional  elastodynamics  problems. 

As  expected,  at  the  same  number  of  dof,  the  increase  in  the  order  of  the  standard  finite 
elements,  the  spectral  elements  and  the  isogeometric  elements  leads  to  the  increase  in  the 
accuracy  of  numerical  results  (however,  the  computational  costs  of  higher-order  elements 
are  greater  than  that  of  lower-order  elements  at  the  same  dof);  see  Figs.  9-10. 

Except  the  known  case  of  the  linear  elements  with  the  lumped  mass  matrix,  other  space- 
discretization  techniques  considered  require  small  time-increments  for  time  integration  at 
the  stage  of  basic  computations.  Moreover,  the  time  increments  should  be  decreased  with 
the  increase  in  the  observation  time.  Even  if  for  explicit  time-integration  methods  the 
time  increments  close  to  the  stability  limit  may  yield  accurate  results  at  a  small 
observation  time,  the  time  increments  should  be  significantly  decreased  with  the  increase 
in  the  observation  time;  see  Fig.  11.  For  example,  for  the  2nd-order  time-integration 
methods,  the  size  of  time  increments  should  be  inversely  proportional  to  the  square  root 
of  the  observation  time  (as  predicted  by  the  exact  time  error  estimator  developed  in  our 
paper);  see  Fig.  1 1 .  Similar  to  explicit  time-integration  methods,  comparably  small  time 
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Comparison  of  accuracy  of  different  order  elements  with  the  consistent 
(non-diagonal)  mass  matrix  that  are  used  for  elastodynamics 
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ID  Impact  Problem:  Velocity  distribution  at  time  T=18. 

Linear  elements  with  reduced  dispersion  (RD)  are  more  computationally  effective  than  other. 

Fig.  9 


Comparison  of  accuracy  of  different  order  elements  with  the  diagonal  mass 
matrix  that  are  used  for  elastodynamics 
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ID  Impact  Problem:  Velocity  distribution  at  time  T=18. 

Linear  elements  with  reduced  dispersion  (RD)  are  more  computationally  effective  than  other. 


Fig.  10 
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The  effect  of  the  size  of  time  increments  in  basic  computations  on  the 
accuracy  of  spectral  10th-order  elements 
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Fig.  11 


increments  should  be  used  with  implicit  time-integration  methods  for  the  1-D  impact 
problem.  According  to  our  results,  the  typical  statement  in  finite  element  textbooks,  that 
for  explicit  methods  a  time  increment  should  be  close  to  the  stability  limit  is  not  true, 
because  the  size  of  a  time  increment  should  depend  on  the  observation  time  and  should 
be  much  smaller  than  the  stability  limit  at  large  observation  times. 


The  comparison  of  the  space-discretization  techniques  based  on  the  non-diagonal  mass 
matrices  (used  with  implicit  time-integration  methods)  show  that  at  the  same  number  of 
dof,  the  isogeometric  elements  yield  more  accurate  results  compared  with  the  standard 
high-order  finite  elements  and  the  linear  elements  with  reduced  dispersion;  see  Fig.  9. 
However,  when  we  compare  the  computational  costs  at  a  given  accuracy,  the  numerical 
results  show  that  the  linear  elements  with  reduced  dispersion  are  more  computationally 
efficient  than  other  space-discretization  techniques.  We  should  also  mention  that 
compared  with  the  isogeometric  3rd-order  elements,  the  computational  efficiency  of  the 
linear  elements  with  reduced  dispersion  decreases  with  the  increase  in  the  observation 
time. 


The  comparison  of  the  space-discretization  techniques  based  on  the  diagonal  mass 
matrices  (used  with  explicit  time-integration  methods)  show  that  at  the  same  number  of 
dof,  the  spectral  high-order  elements  yield  more  accurate  results  compared  with  the 
standard  linear  and  quadratic  finite  elements  and  the  linear  elements  with  reduced 
dispersion;  see  Fig.  10.  However,  when  we  compare  the  computational  costs  at  a  given 
accuracy,  the  numerical  results  show  the  linear  elements  with  reduced  dispersion  are 
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Comparison  of  accuracy  of  different  space-discretization  techniques  used  for 

elastodynamics 

The  size  At  of  time  increments  at  the  filtering  stage  as  an  indicator  of  accuracy 
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(similar  to  the  modal  decomposition) 
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more  computationally  efficient  than  the  spectral  high-order  elements.  We  should  also 
mention  that  compared  with  the  spectral  lOth-order  elements,  the  computational 
efficiency  of  the  linear  elements  with  reduced  dispersion  decreases  with  the  increase  in 
the  observation  time. 


It  is  interesting  to  note  that  the  size  of  time  increments  at  the  filtering  stage  of  the  two- 
stage  time-integration  technique  (this  size  is  calculated  according  the  special  formulas) 
defines  the  range  of  actual  frequencies  used  in  numerical  solutions  and  can  serve  as  a 
quantitative  measure  for  the  comparison  and  the  prediction  of  the  accuracy  and  the 
computational  efficiency  of  different  space-discretization  techniques;  see  Fig.  12. 

The  estimation  and  comparison  of  accuracy  of  different  space-discretization  techniques 
obtained  for  impact  problems  for  which  all  frequencies  of  the  semi-discrete  equations  are 
excited  are  also  valid  for  any  transient  acoustics  or  elastodynamics  problem  for  which 
only  a  part  of  frequencies  of  the  semi-discrete  equations  is  excited.  We  showed  that  the 
two-stage  time-integration  approach  can  be  equally  applied  to  wave  propagation 
problems  under  impact  loading  as  well  as  under  low-  and  high-frequency  loading  (the 
same  range  of  spurious  high-frequencies  should  be  filtered  independent  of  applied 
loading).  The  comparative  study  of  other  space-discretization  techniques  (similar  to  that 
considered  here)  will  help  us  determine  the  most  computationally  efficient  technique  for 
elastodynamics  and  acoustics. 
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Our  publications  with  a  detailed  description  of  the  results  presented  above. 
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4.  Propagation  of  acoustic  waves 


We  have  analyzed  in  more  detailed  the  application  of  the  new  numerical  approach  to  propagation 
of  acoustic  waves.  For  the  space-discretization,  we  have  used  the  standard  linear  finite  elements 
as  well  as  the  linear  elements  with  reduced  dispersion.  However,  any  other  space-discretization 
techniques  such  as  the  spectral  elements,  the  isogeometric  elements,  the  boundary  element 
method,  different  meshless  methods  and  many  others  can  be  applied  with  the  new  approach 
suggested.  The  new  findings  can  be  summarized  as  follows. 

In  contrast  to  the  results  known  from  the  literature,  we  use  the  dispersion  analysis  of  the 
reduced  dispersion  technique  and  the  standard  finite  element  approach  for  transient 
acoustic  problems  in  order  to  detennine  the  effect  of  the  size  of  time  increments  on  the 
dispersion  error.  We  have  found  that  for  linear  finite  elements  the  decrease  in  the  size  of 
time  increments  leads  to  the  increase  in  the  dispersion  error  for  both  techniques. 
Therefore,  the  time  increments  close  to  the  stability  limit  in  basic  computations  yield  the 
most  accurate  numerical  solutions.  We  should  mention  that  these  results  are  totally 
different  from  those  for  multi-dimensional  elastodynamics  problems  for  which  smaller 
time  increments  in  basic  computations  yield  more  accurate  solutions  for  the  linear 
elements  with  reduced  dispersion. 

We  have  shown  that  for  a  suddenly  applied  load  even  the  reduced  dispersion  technique 
with  the  time  increments  close  to  the  stability  limit  leads  to  divergent  results  at  mesh 
refinement;  e.g.,  see  Fig  14  (after  basic  computations).  This  is  explained  by  the  increase 
in  the  amplitudes  of  spurious  oscillations  at  mesh  refinement.  Therefore,  in  order  to  get 
convergent  results,  we  have  modified  and  applied  the  two-stage  time  integration 
technique  to  acoustic  problems.  This  technique  allows  the  quantification  of  the  range  of 
spurious  oscillations  and  their  filtering  from  numerical  solutions  for  different  space- 
discretization  methods  at  the  filtering  stage.  For  the  first  time  we  have  shown  that  the 
application  of  the  two-stage  time  integration  approach  yields  accurate  convergent  results 
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Propagation  of  acoustic  waves  in  the  2-D  case  m 
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[1]  A.  V.  Idesman,  D.  Pham,  Accurate  finite  element  modeling  of  acoustic  waves,  Computer 
Physics  Communications,  185,  2014,  pp.  2034-2045. 


Fig.  13 

Propagation  of  acoustic  waves  in  the  2-D  case 
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The  velocity  distribution  along  AD  at  the  observation  time  T  =  0.8  after  basic 
computations  and  after  the  filtering  stage  using  mesh  refinement  and  the  elements 

with  reduced  dispersion  (MIR). 

Fig.  14 
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at  mesh  refinement;  e.g.,  see  14  (after  filtering).  It  is  also  interesting  to  note  that  along 
with  the  damping  of  spurious  oscillations,  the  filtering  stage  significantly  reduces  the 
numerical  anisotropy  of  solutions;  e.g.,  see  Fig.  13. 


In  paper  by  B.  Yue,  M.N.  Guddati,  J.  Acoust.  Soc.  Am.  118  (2005)  2132-2141,  the 
reduced  dispersion  technique  and  the  standard  technique  with  the  lumped  mass  matrix  are 
compared  at  the  same  time  increments  (corresponding  to  the  Courant  number  x=0.75) 
which  are  close  to  the  stability  limit  for  the  reduced  dispersion  technique  but  are  smaller 
than  the  stability  limit  (x=l)  for  the  standard  approach.  Therefore,  this  comparison  cannot 
be  considered  as  the  illustration  of  the  advantage  of  the  reduced  dispersion  technique 
because  the  time  increments  close  to  the  stability  limit  (x=l)  significantly  improve  the 
accuracy  of  the  standard  approach.  We  have  compared  these  techniques  for  the  time 
increments  close  to  the  stability  limit  for  each  technique.  It  is  interesting  to  note  that 
along  the  Cartesian  axes,  the  reduced  dispersion  technique  and  the  standard  technique 
yield  approximately  the  same  accuracy.  However,  the  numerical  solution  for  the  standard 
technique  is  inaccurate  in  other  directions  as  predicted  by  the  dispersion  analysis.  After 
the  filtering  of  spurious  oscillations  in  all  directions,  we  showed  that  the  reduced 
dispersion  technique  is  much  more  accurate  than  the  standard  approach. 


Our  publications  with  a  detailed  description  of  the  results  presented  above. 

1.  Idesman  A.V.,  Pham  D.  ,  Accurate  finite  element  modeling  of  acoustic  waves,  Computer 
Physics  Communications,  2014,  185,  pp.  2034-2045. 

2.  Idesman  A.V.,  Pham  D.,  Foley  J.,  Schmidt  M.  Accurate  solutions  of  wave  propagation 
problems  under  impact  loading  by  the  standard,  spectral  and  isogeometric  high-order 
finite  elements.  Comparative  study  of  accuracy  of  different  space-discretization 
techniques,  Finite  Elements  in  Analysis  and  Design,  2014,  88,  pp.  67-89. 

3.  Idesman  A.V.,  Pham  D.,  Finite  element  modeling  of  linear  elastodynamics  problems 
with  explicit  time-integration  methods  and  linear  elements  with  the  reduced  dispersion 
error,  Computer  Methods  in  Applied  Mechanics  and  Engineering,  2014,  271  pp.  86-108. 
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5.  Wave  propagation  in  inhomogeneous  materials. 


We  have  extended  the  new  technique  to  wave  propagation  in  inhomogeneous  materials  including 
composite  and  functionally  graded  materials  and  we  have  shown  the  main  differences  for  the 
modeling  of  wave  propagation  in  inhomogeneous  and  homogeneous  materials.  Below  we  will 
summarize  the  main  findings. 


For  the  special  design  of  finite  element  meshes  in  the  1-D  case,  we  can  obtain  very 
accurate  results  for  wave  propagation  problems  in  inhomogeneous  materials  if  we  use  the 
standard  linear  elements  with  the  lumped  mass  matrix  and  the  explicit  central  difference 
method  with  the  time  increments  equal  to  the  stability  limit.  In  this  case,  there  are  no 
spurious  oscillations  in  the  numerical  solutions  after  basic  computations  (the  filtering 
stage  is  not  required)  and  a  discontinuity  in  a  solution  (if  presented)  is  spread  over  three 
nodes.  In  contrast  to  the  finite  element  formulations  with  graded  finite  elements  used  in 
the  literature  for  such  problems,  we  use  a  piecewise  constant  variation  of  materials 
properties  (they  are  constant  within  any  finite  element).  Nevertheless,  even  for 
functionally  graded  materials  the  numerical  results  converge  very  fast  to  the  exact 
solution  at  mesh  refinement.  It  is  interesting  to  mention  that  in  the  1-D  case  the  analytical 
solution  to  wave  propagation  in  functionally  graded  materials  can  be  obtain  by  the 
approach  developed  by  Chiu  T-C,  Erdogan  F.  (One-dimensional  wave  propagation  in  a 
functionally  graded  elastic  medium.  J  Sound  Vib  1999;222(3):453-87).  However,  this 
analytical  solution  is  based  on  the  infinite  series  and  may  lead  to  spurious  oscillations  due 
to  Gibbs  phenomena  when  a  finite  number  of  terms  is  used  in  the  series.  In  contrast  to 
this  case,  the  numerical  solutions  based  on  the  special  non-uniform  meshes  with  standard 
linear  finite  elements,  on  the  lumped  mass  matrix  and  on  the  explicit  central-difference 
method  with  the  time  increments  equal  to  the  stability  limit  do  not  include  spurious 
oscillations  at  any  number  of  degrees  of  freedom.  They  can  be  used  as  reference  solutions 
for  testing  computer  codes  as  well  as  for  the  analysis  of  wave  propagation  phenomena 
under  high-frequency  and  impact  loadings. 

Except  one  case  described  above,  in  all  other  cases  (e.g.,  unifonn  meshes  or  the  non- 
lumped  mass  matrix  or  smaller  time  increments  or  other  time-integration  methods  or 
high-order  finite  elements  or  the  2-D  and  3-D  problems  etc.)  numerical  solutions 
obtained  by  existing  methods  may  include  large  spurious  high-frequency  oscillations 
especially  under  impact  loading.  However,  the  new  numerical  approach  based  on  the 
two-stage  time-integration  technique  (with  basic  computations  and  the  filtering  stage) 
quantifies  and  removes  the  spurious  oscillations  for  these  cases.  For  the  first  time  we 
have  obtained  accurate  finite  element  solutions  of  wave  propagation  problems  in 
composite  and  functionally  graded  materials  without  spurious  oscillations  and  without  the 
interaction  between  user  and  computer  code;  see  Figs.  15-17. 
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Wave  propagation  in  composite  materials  (1-D  case)  1 1 1 

^ _  Non-uniform  meshes  with  the  size  of 
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The  velocity  distribution  along  the  bar  at  the  observation  time  T  =  98.3  (linear  finite  elements). 

(1]  A  V  Idesman,  Accurate  finite  element  modeling  of  wave  propagation  in  composite  and  functionally  graded  materials  Composite  Structures.  117.  2014,  pp.  298-308 

Fig.  15 

Wave  propagation  in  functionally  graded  materials  (1-D  case) 
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V  Idesman,  Accurate  finite  element  modeling  of  wave  propagation  in  composite  and  functionally  graded  materials,  Composite  Structures.  11 7,  2014  pp  298-308 


Fig.  16 
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Wave  propagation  in  composite  materials  (2-D  axisym metric  case) 

Non-uniform  meshes  with  the  size  of 
finite  elements  proportional  to  the  local 
wave  velocity 

^2  Axial  velocity  distribution  along  the  axis  of  revolution  AD  at  T=6  Axial  velocity  distribution  along  the  axis  of  revolution  AD  at  T=6 


Axial  velocity  distribution  along  AD  at  T=6  **\a\  velocity  distribution  along  surface  BC  at  T=6 

Vz  Mesh  refinement  with  non-uniform  meshes  (after  filtering)  Vz  Mesh  refinement  with  non-uniform  meshes  (after  filtering) 


The  axial  velocity  distribution  at  T  =  6  (linear  finite  elements,  averaged  mass  matrix). 

[1]  A.  V.  Idesman,  Accurate  finite  element  modeling  of  wave  propagation  in  composite  and  functionally  graded  materials.  Composite  Structures.  117.  2014,  pp.  298-308. 


Fig.  17 


In  many  known  papers  on  the  simulation  of  the  dynamic  response  of  structures  made  of 
inhomogeneous  materials,  uniform  meshes  or  meshes  independent  of  material  properties 
are  used  with  the  standard  approaches  with  explicit  or  implicit  time-integration  methods. 
We  have  shown  that  despite  unifonn  meshes  are  justified  for  homogeneous  materials, 
non-unifonn  meshes  with  the  constant  ratios  of  the  local  wave  velocity  to  the  size  of  the 
finite  element  yield  much  more  accurate  results  than  uniform  meshes  do;  see  Figs.  15-16. 


We  have  shown  that  for  inhomogeneous  materials  (similar  to  homogeneous  materials), 
the  finite  element  formulation  with  reduced  dispersion  (based  on  the  averaged  mass 
matrix  and  the  linear  elements)  yields  much  more  accurate  results  than  the  standard 
formulations  with  the  lumped  or  consistent  mass  matrices  do;  see  Figs.  15-16.  There  is 
one  exception  for  the  case  of  the  lumped  mass  matrix  with  the  non-unifonn  mesh  with 
linear  elements  and  the  explicit  central  difference  method  with  the  time  increments  equal 
to  the  stability  limit.  However,  this  exception  is  valid  in  the  1  -D  case  only. 


We  have  shown  that  the  size  of  time  increments  at  the  filtering  stage  of  the  two-stage 
time-integration  technique  defines  the  range  of  actual  frequencies  used  in  numerical 
solutions  and  can  serve  as  a  quantitative  measure  for  the  comparison  and  the  prediction  of 
the  accuracy  and  the  computational  efficiency  of  different  space-discretization 
techniques  used  for  wave  propagation  in  inhomogeneous  materials.  The  smaller  size  of 
time  increments  at  the  filtering  stage  calculated  according  to  the  special  formulas 
corresponds  to  more  accurate  numerical  solutions. 
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Our  publications  with  a  detailed  description  of  the  results  presented  above. 

1 .  Idesman  A.  V.  ,  Accurate  finite  element  modeling  of  wave  propagation  in  composite  and 
functionally  graded  materials,  Composite  Structures,  2014,  117,  pp.  298-308. 


6.  Propagation  of  high-frequency  pulses  in  the  Hopkinson  Pressure  Bar 


We  have  considered  the  modeling  of  wave  propagation  in  the  components  of  the  Hopkinson  Bar. 
For  the  first  time,  detailed  accurate  numerical  solutions  for  elastic  wave  propagation  in  a  long 
axisymmetric  elastic  bar  under  impact  loading  are  obtained  using  the  new  finite  element 
technique.  In  contrast  to  known  numerical  techniques,  the  new  numerical  approach  quantifies 
and  removes  spurious  high-frequency  oscillations  which  may  invalidate  numerical  results  in 
impact  loading  simulations.  The  comparison  of  the  accurate  experimental  results  (by  Dr.  Mates 
from  NIST)  for  the  impact  of  striker  and  incident  bars  with  the  corresponding  accurate  numerical 
results  (see  Fig.  18)  allows  us  to  explain  some  details  of  elastic  wave  propagation  in  long  bars. 


Elastic  impact  of  striker  and  incident  bars  (2-D  axisymmetric  formulation) 


4=0.25  m  .  D  4,=  1.5  m  (F 


i 

Axial  strain  evolution  in  the  middle  of  the  incident  bar  at  point  P  (see  [1]) 


[1]  Accurate  finite  element  simulation  and  experimental  study  of  elastic  wave  propagation  in  a  long  cylinder  under  impact 
loading.  Idesman  A.V..  Mates  S..  Impact  Engineering,  2014.  71.  pp.  1-16. 


Fig.  18 


For  example,  due  to  the  absence  of  very  high  frequencies  in  the  obtained  experimental  results, 
the  mathematical  formulation  of  the  problem  should  include  physical  damping  for  the 
corresponding  range  of  high  frequencies.  This  range  can  be  defined  by  the  filtering  stage  of  the 
new  approach  in  terms  of  the  number  of  finite  elements  along  the  radial  direction  of  the  bar.  By 
the  variation  of  this  number  we  can  fit  the  experimental  curves  with  the  numerical  results 
obtained  by  the  new  numerical  technique.  However,  for  the  accurate  numerical  solution  of  the 
impact  problem  with  zero  physical  damping,  the  number  of  elements  in  the  radial  direction 
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should  be  large.  By  the  comparison  of  the  numerical  and  experimental  data,  we  can  accurately 
determine  the  longitudinal  wave  velocity  from  the  experiments.  The  accurate  numerical  solutions 
also  allow  the  analysis  of  the  uniformity  of  the  different  strain  and  velocity  components  across 
the  radius  at  different  distances  from  the  impact  face.  The  validity  of  some  assumptions  used  in 
the  1-D  theory  for  wave  propagation  in  long  bars  has  been  also  checked  by  the  use  of  the 
accurate  numerical  solution.  We  have  also  shown  that  at  the  elastic  impact  the  known  dispersion- 
correction  technique  used  for  the  description  of  the  shape  of  the  wave  pulse  at  different  locations 
along  the  axisymmetric  bar  is  inaccurate  for  the  prediction  of  pulses  close  to  the  impact  face. 

Wave  propagation  in  incident  bar  (2-D  axisymmetric  formulation)  [1] 


ti  t2 

b) 

£z 


a)  c) 

a)  Geometry,  b)  Boundary  conditions  at  the  impact  face, 
c)  Axial  strain  evolution  in  gage  1  for  a  short  time  (just  for  the  first  oscillation). 

[1]  Idesman  A.  V..  Foley  J  .  Dodson  J.  Effect  of  Radial  Confinement  on  Wave  Propagation  and  Vibrational  Response  in  liars. 

Experimental  Mechanics,  2015  (in  preparation). 

Fig.  19 

We  have  also  simulated  the  experiments  on  wave  propagation  in  long  bars  under  impact  loading. 
These  experimental  data  have  been  provided  by  the  group  of  Dr.  Foley  from  AFRL/RWMF, 
Eglin.  The  experimental  results  included  the  evolution  of  the  axial  velocity  at  four  gages  located 
along  an  incident  bar;  see  Figs.  19,  20.  In  order  to  simplify  the  simulations,  we  have  not 
considered  the  contact  interaction  of  a  striker  bar  with  the  long  incident  bar.  The  action  of  the 
striker  bar  we  have  replaced  by  the  following  boundary  conditions  at  the  impact  face:  the  normal 
velocity  vz(r,0,t)  =  Vo(t)  and  zero  tangential  tractive  forces  are  applied  until  time  t  <  t2  (we 
matched  the  numerical  and  experimental  results  for  the  first  oscillation  in  gage  1  by  an 
appropriate  selection  of  Vo(t);  see  Fig.  19)  and  the  zero  normal  and  tangential  tractive  forces  are 
applied  after  time  t  >  t2.  The  problem  was  solved  by  the  new  numerical  technique.  We  have 
obtained  a  surprisingly  good  agreement  between  the  experiments  and  the  calculations  not  only 
for  large  oscillations  but  also  for  the  following  small  oscillations;  see  color  and  black  curves  in 
Fig.  20.  We  should  emphasize  that  in  our  approach  we  assigned  the  velocity  Vo(t)  at  the  impact 
face  that  allowed  us  to  fit  the  experimental  and  numerical  results  just  for  the  first  big  oscillations 
in  gage  1  for  a  short  time;  see  Fig.  19.  Then,  all  other  numerical  results  at  time  t  >  0.00015  s  in 
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gage  1  and  the  results  in  gages  2-4  represent  an  accurate  numerical  solution  of  the  wave 
propagation  problem  for  the  given  boundary  conditions;  see  Fig.  20.  By  the  numerical  modeling 
we  could  explain  such  interesting  experimental  results  as  a  small  increase  in  the  amplitudes  of 
big  oscillations  as  well  as  the  appearance  of  new  small  oscillations  with  the  increase  in  the 
observation  time  (the  effect  of  high  frequencies  and  physical  dispersion);  see  Fig.  20. 


Wave  propagation  in  incident  bar  (2-D  axisymmetric  formulation)  [1] 


The  axial  strain  evolution  in  gages  1  (a),  2  (b),  3  (c)  and  4  (d)  for  a  solid  bar.  Black  curve  -  numerical  simulation, 

color  curves  -  experiments. 

[  1]  Idesman  A.V.,  Foley  J..  Dodson  J.  Effect  of  Radial  Confinement  on  Wave  Propagation  and  Vibrational  Response  in  Bars. 

Experimental  Mechanics,  2015  (in  preparation). 

Fig.  20 


Our  publications  with  a  detailed  description  of  the  results  presented  above. 


1.  Idesman  A.V.,  Mates  S.  P.,  Accurate  finite  element  simulation  and  experimental  study  of 
elastic  wave  propagation  in  a  long  cylinder  under  impact  loading,  International  Journal  of 
Impact  Engineering,  2014,  71,  pp.  1-16. 

2.  Idesman  A.V.,  Foley  J.,  Dodson  J.  Effect  of  Radial  Confinement  on  Wave  Propagation 
and  Vibrational  Response  in  Bars,  Experimental  Mechanics,  2015  (in  preparation). 
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7.  Activities  supported  by  the  grant 


Students 
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